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Abstract. Numerical simulations are performed of the gravitational collapse of 
a scalar field with a \<j) A potential. Comparisons are made with the thin shell 
approximation. 
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1. Introduction 



Domain walls in general relativitypQ are usually treated in the thin shell 
approximation. [2] However, there is also a description of a domain wall as a soliton 
of a field theory. In a certain limit one expects the field theory description to reduce to 
the thin shell description. Indeed, one can formally expand the Einstein-scalar equations 
in powers of the thickness and obtain the thin shell equations at lowest order in this 
expansion. [3111] Nonetheless, one cannot simply assume the validity of such an expansion. 
Instead in a field theory one gets to choose initial data and then the field equations give 
the results of the evolution of those data. Thus the thin shell approximation is a good 
approximation for a thick domain wall to the extent that the evolution of initial data 
that are well approximated by the thin shell treatment continue under evolution to be 
well approximated by this treatment. To test this approximation, we will choose initial 
data for a spherically symmetric thick domain wall, perform a numerical evolution, 
and compare to the corresponding thin shell approximation. Methods are described in 
section 2, results in section 3 and conclusions in section 4. 



2. Methods 



The Lagrangian for a thick wall takes the form 

L = -|V a 0V a 0-l/(0) (1) 

where V{<p) is a potential with two minima. The equation of motion associated with 
this Lagrangian is 

V a V a - ^ = (2) 
The shell also satisfies the Einstein field equation 

G a b = KT a b (3) 

Here G a b is the Einstein tensor, k = 8irG where G is Newton's gravitational constant, 
and the stress-energy of the scalar field is given by 

Tab = V a 0V b - |<u(V c 0V c + 2V) (4) 

Simulations of such thick walls in the case where their self-gravity can be neglected were 
done by(5] and [Bj. We will use the standard A0 4 potential 

V = \tf-r?f (5) 

For simulations in spherical symmetry, one often takes the area radius as one of 
the spatial coordinates and chooses time to be orthogonal to this radial coordinate. [7] 
However, this coordinate system breaks down when a trapped surface forms and thus 
cannot follow the evolution after black hole formation. Instead we will use the method of 
[8] and use maximal slicing with radial length as the radial coordinate. Maximal slicing 
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allows us to simulate part of the region inside the black hole without encountering the 
singularity. The metric takes the form 

ds 2 = -a 2 alt 2 + {dr + r dtf + R 2 {d6 2 + sin 2 #V) (6) 

Note that the usual area radius R is not one of the coordinates and is instead a function 
of the coordinates t and r. Thus the spatial metric i a b has components 

Irr = 1 

lee = R 2 

7w = R\m 2 9 (7) 

The extrinsic curvature, K ab is defined by 

K ab = -ia C V c n b (8) 

where n a is the unit normal to the surfaces of constant time t. However, due to spherical 
symmetry and maximal slicing, there is only one independent component of the extrinsic 
curvature. Specifically we have 

K\ = K% = -\K\ (9) 

Equation (El) is equivalent to 

d tlij = -laKij + + Djfc (10) 

where Di is the covariant derivative of the spatial metric jij. The rr component of eqn 
(HUD yields 

d r (3 r = aK r r (11) 
whose solution is 

P r = [ r aK r r dr (12) 
Jo 

The 99 component of eqn (TTUT) yields 

OL 

d t R = (3 r d r R + -RK r r (13) 

We now use the momentum constraint of the Einstein field equation to determine 
the extrinsic curvature. For maximal slicing (K = 0) this constraint is 

D a K ab = -Ki bc n d T cd (14) 
Define the quantities P and S by 

P = n a V a 0, S = d.,4 (15) 
Then eqn (T141 becomes 

d r K r r + ZR^K^r = -kPS (16) 

Note that there is also a Hamiltonian constraint associated with the Einstein field 
equation. In the case of maximal slicing, this constraint is 

^R-K ab K ab = 2KT ab n a n b (17) 
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where is the spatial scalar curvature. This equation yields 



d r d r R 



1 - (d r R)' 
2R 



\R\MK r r f + K{P 2 + S 2 + 2V) 



It 



We now determine the lapse a. It follows from the maximal slicing condition that 



D a D a a 



a 



K ab K ab 1 ''" ;< 



l-^(«V + 7 » 6 ) 



which yields 



d r d r a + —(d r R)(d r a) 
R 



a: 



(K%) 2 + k{P 2 - V) 



(19) 



(20) 



We now consider the evolution of the scalar field. From the definitions of P and S 
it follows that 

d t (j) = aP + f3 r S (21) 
d t S = a(d r P + K T r S) + Pd r a + (3 r d r S (22) 

The equation of motion, eqn (j2J) becomes after some straightforward but tedious algebra 



8 t P = p r d r P + Sd r a + a 



R- 2 d r (R 2 S) - %■ 



(23) 



The initial data are as follows: we choose a moment of time symmetry so that P 
and K a b vanish. The scalar field is chosen to have a flat spacetime domain wall profile 
with the center of the wall at a radius r . More precisely, define the quantities e and a 
by 

1 



a 



(24) 

lV2Xr] 3 (25) 

where rj and A are the parameters of the potential. Here a is the energy per unit area of 
the wall and e is an effective wall thickness. Formally the thin shell limit of the solution 
is the limit as e — > at constant a. Correspondingly, one can specify e and a; then A 
and r) are determined by 
2 



A 



3ae 3 



rere 



(26) 
(27) 



In flat spacetime a static planar domain wall solution is given by 

= ?7tanh(z/e) (28) 
We choose for initial data 

= ?7tanh((r — r )/e) (29) 

where the constant r$ can be chosen arbitrarily but should be chosen to be much larger 
than e. The quantity S is set equal to <9 r 0. Eqn. (fT8"j) is integrated for R using the fact 
that at the origin R = and d r R = 1. Then at each time step, the evolution proceeds 
as follows: first eqns ( TL61) and ( fTBl are integrated to find K r r and R. And then eqn (1201) 
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Figure 1. The scalar field <j> at times 0, 2, 4, 6 and 



is solved for the lapse a (using a tridiagonal method and the fact that d r a = at the 
origin and a — > 1 at infinity). Then eqn ( fl2|) is integrated to find the shift /3 r . Finally, 
the quantities (ft, S and P are evolved to the next time step using eqns ( l2il . ( )22l) and 
( |23|) respectively. The evolution is done using the iterated Crank-Nicholson method, 
and all spatial derivatives are found using standard centered differences. We use units 
where k = 1. 

As the evolution proceeds we can check for black hole formation by looking for the 
presence of a marginally outer trapped surface. In spherical symmetry, such a surface 
is given by the condition V a -RV a -R = which in our coordinate system is equivalent to 

d r R + \RK% = (30) 



3. Results 



The simulations are run with n + 1 points evenly spaced between r = and a maximum 
value r max . We choose r max = 20. We would like to know how the initial scalar 
field profile changes under evolution. Fig [1] shows the result of a simulation with 
r = 5, e = 0.25 and a = 0.15. This simulation was run with n = 9600 and results 
are shown for the scalar field <fi at times 0, 2, 4, 6, and 8. As the evolution proceeds, 
the scalar field profile moves inward. Furthermore, the profile becomes steeper, and 
departures from the simple tanh form become more pronounced. A marginally outer 
trapped surface forms at t = 5.09. 

In order to be sure that these results are not numerical artifacts, we need to know 
that the code is convergent and that these results are within the convergent regime. 
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Figure 2. ln||C|| vs. t for n — 9600 (solid line) and n 
results demonstrate second order convergence. 



19200 (dashed line). The 



Note that eqn. (|T3|) is not used in the evolution, but should nonetheless be satisfied to 
within numerical error due to finite differenceing. Thus this equation provides a check 
on the performance of the simulation. More precisely, define the constraint quantity C 
by 

C = ^K r r + R~ 1 {f3 r d r R-d t R) (31) 

and let ||C|| be the L 2 norm of C. Then in an exact solution C should vanish, while in a 
numerical treatment C should converge to zero in the limit of zero step size. Fig [2] shows 
ln||C|| as a function of time. Here the parameters are as in the previous simulation, 
except that one simulation is run with n = 9600 and one with n = 19200. The results 
demonstrate second order convergence: that is, halving the step size reduces C by a 
factor of 4. 

In the thin shell formalism one is mostly concerned with the motion of the wall. 
The motion of a spherical wall is described by giving its area radius as a function of 
proper time. That is, 

R = Ro(r) (32) 

Using the results of [TJ [2] one can show that the equation of motion of the thin shell is 

R = |k<t(1 + Rlf 2 - 2^(1 + Rl) (33) 

We would like to know how well the motion given by solving eqn. (|33|) models the 
behavior of the scalar field domain wall. Here we can do a direct comparison: at any 
given time, the position of the wall will be taken to be the place where = and the 
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proper time r will be that of an observer who is always at the position of the wall. At 
each time step of the simulation, we can find the position of the wall, so it remains to 
evaluate r. Let u a be the four-velocity of the observer who remains at the position of 
the wall. Then it follows that u a is a unit timelike vector for which u a V a <P = 0. It then 
follows using eqns (JHJ) and ( TT5l) that 



u 



A' 



n 



Pfd_ 
S \dr 

where the quantity A is defined by 



.4 




(34) 



(35) 



evaluated at the position of the wall. The relation between the normal vector and the 
time coordinate is 



n a = -aV a t 
so using eqn. (CM|) we have 

dt „„ 
-r = u a V a t -- 
ar 

and therefore we find 

dr = aAdt 



-a 1 u a n a 



a 



-i 



A' 



(36) 



(37) 



(38) 



Choosing r = at the begining of the simulation, we then integrate eqn. ( 138]) to find r 
at all times of the simulation. Since we also have Rq at all points of the simulation, we 
can produce the thick wall Rq(t) for comparison with the Rq(t) of the thin wall. Fig. 
([3]) shows such a comparison. Here the parameters of the wall are the same as in the 
previous simulations, with the solid line being that of the simulation and the dashed 
line the solution of eqn. ( |33l) . 

We now consider a comparison of the stress-energy of the simulation to that of 
the thin shell approximation. From eqns (jl]) and (fT5|) it follows that the trace of the 
stress-energy tensor is 

T = P 2 — S 2 — AV (39) 

However, in the treatment of [3] to lowest order in thickness of the wall we have 

-9a 



T 



(40) 



4ecosh 4 (;z/e) 

where z is geodesic distance from the center of the wall. Calculating geodesic distance 
from a simulation can be complicated; but luckily because the expression for T in eqn. 
( 1401) falls off so rapidly, we only need z in the vicinity of the wall's center and thus can 
approximate z using a Taylor series. In particular, note that at the center of the wall 
V a z must be a unit vector orthogonal to u a . It then follows that 

8 V P 



W a z = A~ 



dr 



S 



n 



(41) 
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Figure 3. Rq vs. r for the simulation (solid line) and the thin wall equations of 
motion (dashed line). 



A- 1 



(42) 



Taking the inner product of eqn. (141 p with (d/dr) a we find that at the center of the 
wall 

dz 
dr 

and we therefore find that near the wall z is well approximated by 

z = (r - Ro)/A (43) 

Thus the thin wall expression for T is given by eqn. (141) I) with z given by eqn. (1431 . 

Such a comparison is given in Fig (j4J) for t = 2 and in Fig ([5]) for t = 6. The 
parameters of the wall are the same as in the previous simulations. Here the trace of 
the stress-energy T is plotted for the simulation (solid line representing the expression of 
eqn. (1391) ) and for the thin shell approximation (dashed line representing the expression 
of eqn. (j40l) . Note that the simulation results are well approximated by the thin shell 
expression, though the approximation becomes less good as the evolution proceeds. 



4. Conclusions 



Our simulations indicate that the thin wall approximation is an excellent approximation 
for thick wall gravitational collapse. Perhaps more importantly, we have developed a 
robust numerical method for thick wall collapse that could be used on other projects. 
In particular, we plan to simulate the collapse of a charged thick domain wall. Thin 
charged walls collapse to charged black holes and can even be used to form an extreme 
(charge equals mass) black hole. [9] In contrast, simulations of the collapse of a free 



Gravitational Collapse of Thick Domain Walls 9 
0.2 I 1 1 1 1 1 1 1 1 



-0.2 - 
-0.4 - 




. 1i4 I 1 1 1 1 1 1 1 1 

3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 

r 

Figure 4. T vs. r at t — 2 for the simulation (solid line) and the thin wall 
approximation (dashed line). 
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Figure 5. T vs. r at t — 6 for the simulation (solid line) and the thin wall 
approximation (dashed line). 
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charged scalar field yield black holes that are always far from extreme. [TD] It will be 
interesting to see whether an extreme black hole can be formed by the collapse of a 
charged thick domain wall. 
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